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Abstract 

We compare the grand canonical partition function at fixed chemical potential a 
with the canonical partition function at fixed baryon number B, formally and by 
numerical simulations at u = and B = with four flavours of staggered quarks. 
We verify that the free energy densities are equal in the thermodynamic limit, and 
show that they can be well described by the hadron resonance gas at T < T c and 
by the free fermion gas at T > T c . 

Small differences between the two ensembles, for thermodynamic observables char- 
acterising the deconfinement phase transition, vanish with increasing lattice size. 
These differences are solely caused by contributions of non-zero baryon density sec- 
tors, which are exponentially suppressed with increasing volume. The Polyakov loop 
shows a different behaviour: for all temperatures and volumes, its expectation value 
is exactly zero in the canonical formulation, whereas it is always non-zero in the 
commonly used grand-canonical formulation. We clarify this paradoxical difference, 
and show that the non-vanishing Polyakov loop expectation value is due to contri- 
butions of non-zero triality states, which are not physical, because they give zero 
contribution to the partition function. 
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1 Introduction 



To simulate QCD thermodynamics on the lattice £Q, one commonly uses the grand 
canonical partition function with respect to the quark number as a function of a 
chemical potential \i. It is given, after integration over the fermion fields, by 

Z GC (T, /*) = J [DU] e~ s ^ det M(U; /*) . (1) 

Recently, another approach using a canonical formalism has been used The 
canonical partition function, which will be derived in detail in the next section, is 

Z C (T, Q) = T d (!£) e-^Z GC (T, /x = i^) . (2) 

The physics described by both ensembles, grand-canonical and canonical, must be 
the same in the thermodynamic limit, ie. the free energy density should be the same. 
This has been shown in Ref. jU Ej, and will be confirmed in this study. It is thus 
puzzling that the Polyakov loop expectation value in the canonical ensemble is ex- 
actly zero, while it is non-zero in the grand-canonical ensemble, for all temperatures 
and volumes. We will show that this discrepancy is due to contributions from the 
canonical sectors with a quark number that is not a multiple of three: the so-called 
non-zero triality sectors 3 . Discussions on the role of these non-zero triality sectors 
have a long history [Oj and include speculations that their influence persists even in 
the thermodynamic limit, so that they must be explicitly projected out. It is our 
purpose to clarify these issues. 

In the following, we discuss properties of the grand canonical partition function as 
a function of an imaginary chemical potential and construct the canonical partition 
function (Section II). We then show that the Polyakov loop vanishes in the canonical 
ensemble (Section III), and resolve the paradox above (Section IV). After presenting 
our numerical method to simulate the zero baryon density sector, which is the first 
step towards finite density QCD simulations (Section V), we elaborate on the results 
for the free energy density as a function of the imaginary chemical potential, which 
we compare with predictions of the hadron resonance gas model [7j and of a free 
fermion gas (Section VI). We further study the expectation values and finite size 
effects of thermodynamic observables, like the plaquette and the chiral condensate 
in both formulations (Section VII). Conclusions follow. Preliminary results of this 
study have been presented in Ref. 

3 Triality is defined as the difference between the number of quarks and the number of anti- 
quarks modulo 3. 
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2 Canonical Ensemble 



Let us first discuss symmetries of the grand canonical partition function Z GG (T,fi) 
as a function of an imaginary chemical potential /i = ifii, following Ref. [S|: 

• Evenness: Z GC (-ifii) = Z GG (+i/j,i). 

The transformation /i — > —fi can be compensated by time-reversal, ie. by 
interchanging particles and anti-particles. Time reversal is equivalent to CP 
symmetry (since CPT is always a good symmetry), and thus does not change 
the thermodynamics in the absence of CP violating terms. 

• ^-periodicity in /i 7 : Z GC {i{^i + ^)) = Z GC {ini). 

A shift in the imaginary chemical potential [xi — > \ij + can be exactly 
compensated by a ^-transformation, ie. by a rotation of the Polyakov loop 

■ 27rfc 

Pol(x) — > e~ l ~ Pol(x), configuration by configuration. Since the path integral 
sums over all possible gauge fields, the partition function stays the same. 

These properties lead to an expectation for the phase structure in the T — \xi plane, 
see Ref. jH] and Fig. ^ In the low-temperature phase where the ^-symmetry is 
realised ("disordered"), the periodicity in fij is smoothly realised. However, at high 
temperature, this ^-symmetry is spontaneously broken in favor of one of the Z3- 
sectors ("ordered"). Which sector is preferred depends on fij. Thus, the existence 
of ^-sectors forces the appearance of phase transitions at high temperature. By 
symmetry, these discontinuities must appear at /i/ = ^sML + The transitions are 
first-order since the derivative of the free energy with respect to \i\ is discontinuous. 

To obtain the canonical partition function Z G (T,Q), we fix to Q the conserved 
charge N = J d 3 x ip{x) 70 ip{x), which represents the net quark number. This is 
accomplished by inserting a 5-function in the grand canonical partition function: 



Z C (T,Q) = / [DU][D*][DV] e -s 9 lU;T]-s F im*;T] s (ft-Q\ . 



(3) 



The ^-function admits a Fourier representation with the new variable /^/: 




x e 



S g [U;l3\-S F [U^M+Vi I d 3 xj(x) 70 ip(x) 




x e 



S g lU;f3]-S F [U^M+ifiiT / T drfd 3 x^(x) 7o i>{x) 



(4) 
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Figure 1: Phase diagram of Zcc^H-i) m the (/x/,T) plane. The arrows indicate 
the orientation of the Polyakov loop. The vertical lines mark the "order-order" Z 3 
transitions, which are first-order. Properties of the "order-disorder" Z 3 transitions 
(curved lines) depend on the parameters (number of flavours, quark masses) of the 
theory. 



where Af is a normalisation factor. In the last line, we have used the fact that Q is 
conserved. 

One recognises ifii = iftjT as an imaginary chemical potential, so that 



Z C (T, Q) = M / dfije-^ZadT, ifrT) 



2tt 



dflje-^ZGci^ifljT) , (5) 



where we have exploited the ^^-periodicity in fij of Z GC {i\ii) in the last step 4 . 
From this periodicity, it follows that Zc(T, Q) = except for j = B G Z, where 
I? is the baryon number. The canonical partition function vanishes for non-integer 
baryon number, i.e. for non-zero triality sectors. Note that our argument does 
not rely on particular spatial boundary conditions. The same conclusion holds for 
periodic or, for example, C-periodic spatial b.c, even though the latter break the 



4 Note that the evenness of Zcc{ini) in implies Zc(T,Q) = Zq(T, —Q). In particular the 
Zc(T, Q)'s are real as expected. 
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Z(3) symmetry. For convenience, we write 



ZGc{T,ini) . 



(6) 



Now, the fugacity expansion allows us to go back from the canonical partition 
functions to the grand canonical one: 



This expression is identical to Eq.((TJ), as can be seen by substituting © for Zc(T, Q) 
above and summing over Q first. However, we can remove from the sum the non- 
zero triality (fractional B) sectors, since each of them gives a zero contribution, thus 
obtaining: 



Let us stress again that this grand canonical partition function is identical with 
the one given in Eq.Q]). The two expressions differ by the inclusion of non-zero 
triality sectors, which we just showed are zero. However, these zero contributions 
are explicitly projected out in (jSJ) [U]. In the following, we will come back to this 
and consider the effect of this projecting out on the Z 3 -sensitive Polyakov loop. 

3 Polyakov Loop in the Canonical Ensemble 

The expectation value of the Polyakov loop in the canonical ensemble is zero for 
all temperatures and volumes. We show this explicitly as follows. The chemical 
potential is introduced on the lattice as an external imaginary gauge field 




(7) 



Q=—oo 





B=-oo 




(9) 
(10) 



or equivalently as 





(12) 



on a given temporal hyperplane x^ Q . An imaginary chemical potential i/ii 
can then be absorbed in a Z 3 centre transformation 



• 2nTk 
C 3 



U±{XiX4 = x 4o ) — > e ta 3 Ui(x) = zih^UijyX.x^ = x io ) . 



(13) 



4 



with z(k) = e l ~ I3. As a consequence, the two configurations {[/,///} and 
{z(k)U,fj,i — ^p} have the same value for the Dirac determinant det M(U; fx 1) = 
det M(z(k)U ; fxi — ^f-^) , but the Polyakov loop is centre-rotated. We can then group 
the configurations of a canonical ensemble in triplets having ^-rotated Polyakov 
loop: 



Zc(T,B) = - 



e -i3B^(- 



1 z 

[DU] e - Sa[u '® -^det M(z(k)U 4 (x A = X4o),/x/) 
3 



fc=0 



(14) 

The three members of a triplet give identical contributions to Zc{T,B), since 
det M(z(k)U; m) = det M(U; fxi + ^p) and e"* 3 ^^ = 1 for B e Z. In Fig. 
we show the distribution of the Polyakov loop, where the Z3 symmetry (and hence 
the triplets) is clearly visible in the canonical ensemble (bottom). In each triplet, 

the average of the Polyakov loops is Poli x ^1 + e~ % ^~ + e l ir ) = 0, and therefore 
the ensemble average also vanishes: 



{Pol)z c (T,B) = 



(15) 



for any integer baryon number and temperature. Note again that the argument does 
not depend on a particular choice of spatial boundary conditions. 



4 X 4, T < T r 



A X 4, T > T r _ 



6 X 4j T < T c 



6 X 4, T > T r 



Grand canonical 
ensemble 



Canonical 
ensemble 




Figure 2: Distribution of the complex Polyakov loop trace in the grand canonical 
(top) and canonical (bottom) ensembles, left: 4 3 x 4, right: 6 3 x 4. In the ther- 
modynamic limit, the distributions agree for both ensembles, up to two additional 
^-rotations in the canonical ensemble. 
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4 Polyakov Loop in the Grand Canonical Ensem- 
ble 



In the ensemble generated by the grand canonical partition function Eq.(fT|) 

Z GC (T, fi) = J [DU}e~ s ^ det M(U; /x) , (16) 

the fermion determinant explicitly breaks the Z3 symmetry, so that the expectation 
value of the Polyakov loop 

(PoI)z gc{ t,») ^ (17) 

is non-zero for any chemical potential, temperature and volume. In the following, 
we show that this non-vanishing value is caused by canonical sectors with quark 
numbers which are not a multiple of three. 

We express the grand canonical partition function via the fugacity expansion in 
the quark number Q (we take \x = for notational simplicity only; the argument 
holds for any /i): 

Z GC {T, fi = 0) = J2 Zc{T, Q) with Z C (T, Q) = if Q £ mod 3 
Q 

= ... + Z c (0) + 2c(l) + % c (2) + Z c {3) + Z C (A) + ... , (18) 

where Zc(-) indicates Zc(-) = 0. The canonical partition functions can be written 
as Zc(T,Q) = J2iWi(Q), where i labels each configuration, and Wi{Q) is the cor- 
responding Boltzmann weight. The expectation value of the Polyakov loop is then 
generically given by 

/p „ Eg Num(Q) 

<P0 ' >GC = Z GC (T,, = ) < 19 > 
... + Num(0) + Num(l) + Num(2) + Num(3) + Num(4) + .. 
= ... + Z c (0) + %c(l) + ^c(2) + Z c (3) + J c (4) + ... + ' 

where Num(Q) = J2i P°hWi(Q), which vanishes if Q is a multiple of 3 due to 
Eq.([15|). It follows that the contributions of canonical sectors with fractional baryon 
number to the Polyakov loop are unphysical, since the corresponding canonical ex- 
pectation value is infinite: 

. . Num(Q 7^ mod 3) non-zero , . 

(Pol)z ci T, Q ^ m od 3) = Mc{Q ^ 0mod3) = = °° • (2°) 

This argument makes it clear that the expectation value of the Polyakov loop is 
non-zero for all temperatures and volumes, if we use the grand canonical partition 
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function with respect to the quark number, see Eq.flQ and Fig.|21 top row. However, 
if we use the grand canonical partition function with respect to the baryon number, 
see Eq.fjHJ), the expectation value of the Polyakov loop will be exactly zero even in 
this equivalent grand canonical formulation. 

Thus, the physical meaning of the Z 3 -sensitive Polyakov loop expectation value 
is rather limited. It is the Z 3 -invariant correlator (Pol(0)Pol(xy) which is physical, 
and indicates confinement or deconfinement by its |x| — > oo limit. In the canonical 
ensemble, this limit is not equal to \ (Pol)\ 2 , which is identically zero: the clustering 
property is not satisfied. This is evidence of spontaneous breaking of the center 
symmetry 5 in the presence of fermions. The symmetry is broken spontaneously at 
all temperatures and densities, rather than explicitly as in the usual grand canonical 
ensemble. 



5 Numerical Approach to Zero Baryon Density 

In order to design an algorithm that is able to measure an observable as a function 
of the quark, or rather baryon number, we need to understand how the expectation 
value of an observable O can be evaluated in the canonical ensemble. It is given by 



, as h £r d Pi e ~ t3Bfl1 f[DV) e ~ S9[UM det M (^; m = mT) 6{u) 

{0)b = wj • (21) 

We recognise the following numerical description. We treat dynamical degree 

of freedom, and supplement the ordinary Monte Carlo (Hybrid MC, R-algorithm, 
PHMC, RHMC, . . .) at fixed Jii with a noisy Metropolis update of /Zj — > Ji'j keeping 
the configuration {U} fixed. Thus, we alternate two kinds of Metropolis steps. 

(i) Update of the links by standard Hybrid Monte Carlo : 

Keeping the imaginary chemical potential \ij fixed, we propose a new config- 
uration {U'}, obtained by leapfrog integration of Hamilton's equations, as a 
Metropolis candidate. It is then accepted with the ordinary Metropolis prob- 
ability 

Prob(U -> U') = min (l, e" A5 ) , (22) 
where AS is the difference between the action of {U'} and that of {£/}. 

(ii) Metropolis update of the imaginary chemical potential by a noisy estimator: 
Keeping the gauge field configuration {U} fixed, we propose a new imaginary 



3 We are grateful to L. Yaffc for pointing this out to us. 



7 



chemical potential fi'j, obtained from fij by a random step drawn from an even 
distribution. The update is based on the acceptance 



e -i3B# det Nf (Ipi^'j) + m 
e -i3B Pl det*' (#(///) +m) 



Prob (w - A ) = min ( 1, ., m " Ir T ' I ■ (23) 



The ratio of determinants is evaluated with a stochastic estimator (see Ap- 
pendix EJ), namely 

det^Q^ + m) 1 1 J 

where 77 is a Gaussian complex vector. Since one Gaussian vector is sufficient, 
the computational overhead is negligible. 

The algorithm above allows the sampling of any positive measure in /xj. However, 
the oscillatory part e _ * 3B w j n fa e sampling weight causes a sign problem for non-zero 
baryon number B. One can use | cos(3.B/Xf)| as a sampling measure and fold the 
remaining sign in the observable. But such an approach breaks down for rather small 
B already ^H]- Nevertheless, at B — 0, e _i3S ^ 7 |s=o — 1 an d thus, the Boltzmann 
weight is real and positive. Thus, this simple algorithm suits our purpose here. To 
help ergodicity, we can also perform a "Z3-move" at any time: 

/xj— >/Xj± — — U A {x, x 4 = Xi ) -> U 4 (x, x 4 = x^e^'f, Vx , (25) 

where UA the temporal links at a given time-slice X4 . Such a "Z3- 

move" is always accepted, since the configuration {U, /x/} and the one with a centre- 
rotated Polyakov loop, but shifted imaginary chemical potential, {U x e~ % ~ , /xj + 
^y-} have the same Dirac determinant, and thus the same sampling weight, as dis- 
cussed at length in section El 

A computational detail: For T > T c , the /^-distribution is sharply peaked around 
0,±^p. To sample this distribution accurately in the whole interval, we apply a 
multicanonical algorithm in the T > T c regime for the larger lattices (6 3 x 4 and 
8 3 x 4) jTTj. For this, we bias the sampling of the imaginary chemical potential by 
modifying the acceptance probability 

Prob multi (/x/ -> 6) = min ( 1, det " 7 (ff 04) +™) ^Uas^-Uas^ (26) 
^ W V det^(^(/xz) + m) ; 1 J 
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with bias(ni) chosen such that the sampled histogram becomes flat for all /i/ 6 . The 
expectation value of an observable O is then given by 

(O) = y ^^7J £ d ^ ^~ Uas{Uw) . (27) 

where {£/;/if} labels the configurations {U} sampled at imaginary chemical poten- 
tial fij. 

We focus on four flavours of Kogut-Susskind fermions with degenerate mass 
am = 0.05 and lattices with N t = 4 time-slices, i.e. ¥p = 0.2. With these parameters, 
the zero-temperature pion mass is about 350 MeV ^2]- Simulations are performed 
on lattices with spatial extents 4 3 , 6 3 and 8 3 at seven temperatures 7 , ranging from 

7jr = 0.85 to 1.1, with good overlap between the "neighbouring" ensembles. We 

c 

analyse the results using Ferrenberg-Swendsen reweighting [13] . 



6 Free Energy Density 

In the grand canonical ensemble, the change of the free energy density with chemical 
potential \x (as a dimensionless quantity) is given in terms of the partition function 

AF(T,») _ 1 , Zgo(T,ri 

VT 4 VT* S Z GC {T,0) ' 1 } 

A standard approach [TH ITo] is to perform a Taylor expansion in \x about \i — 0, 
where the derivatives entering the series may be expressed as complicated expecta- 
tion values evaluated at \i — 0. Remember that this expansion is in even powers of 
/i, since Zcc(n) = Z GC (—^). In our approach, the free energy density comes for free 
from the /^/-histogram in the canonical simulation, and moreover, to all orders. At 
low temperature, however, the histograms are quite noisy. Therefore we will, when 
needed, use results from the more sophisticated method [H] , where we can calculate 
the grand canonical partition function for an arbitrary imaginary chemical potential 
as a consequence of the reweighting method that we apply. 

In Fig. 01 we show the free energy divided by VT A versus ^ for |r < 1, ^ ~ 1 
and if- > 1. In all cases, we observe a minimum at ^ = 0. Therefore, in the ther- 
modynamic limit, only = mod -y will survive. This establishes numerically 

6 A simple way to get an estimate of the function bias(fi]) is the following: One starts by 
sampling with no bias to produce a histogram hist(fij) of the sampled /i/. One then fits bias((j,i) 
to — \og(hist(fii)) with a suitable Ansatz like afij — b[ij, or uses a table. 

7 We relate the coupling j3 to the temperature T via T = a ^ Nt an d the perturbative two-loop 
/^-function. 
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T7T C ~0.9 T/T c ~1.0 T/T c -1.1 




■j/3 -it/6 it/6 it/3 -it/3 -it/6 it/6 it/3 -n / 3 -it/ 6 it/6 it/3 



H,/T H,/T n,/T 

Figure 3: Ai ^f ^ as a function of ^f, at temperatures ~ 0.9, 1.0, 1.1 from left to 
right. The free energy density varies much more upon entering the high-temperature 
phase, and the Z 3 first-order transitions become visible (right). 
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Figure 4: AF ^^^ as a function of for jr ~ 0.9. The histogram method is very 
noisy. We show (left) a rescaled version of the leftmost plot in Fig. El We also present 
(right) results based on a reweighting method with variance reduction [3]. The 
results are in agreement with the histogram method, but allow for a more reliable 
description by a Fourier expansion. One Fourier coefficient suffices to describe the 
data points. The reweighting method j3] calculation is computationally demanding 
and has not been performed yet for the 8 3 x 4 lattice. We thus only draw the fit, 
which is based on histogram data. 



the expected equivalence of Zc(T, B = 0) with Zgc(T,[i = 0). 

For y~ ~ 0.9, no singularities develop at ^ = ±| in the thermodynamic limit, 
thus indicating a crossover, as expected from the phase diagram T-fij, Fig. ^ In 
Fig.HJ (left), we show the free energy density, determined by the histogram method, 
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which is very flat and noisy unfortunately. The periodicity of the free energy density 
is ^p, and we exploit it by a Fourier expansion in 3fc^- using the Ansatz 

^fi) = c (l-co S (3^)) + dcos(6^) + ... . (29) 

In order to improve the determination of the coefficients c,d,. . ., we use results based 
on the reweighting method described in Ref. [3]. Within errors, the free energy den- 
sity is in agreement with the histogram method, but with much smaller statistical 
uncertainty. The fit is excellent already with one Fourier coefficient, with no indi- 
cation for higher Fourier components, at least on the small lattices we consider. 

In the hadron resonance gas model (see Ref. ^7] for a detailed discussion), the 
partition function can be split into mesonic and baryonic contributions. Since we are 
interested in the dependency on a baryon chemical potential /i# = 3/i, it is sufficient 
to study the baryonic part only. In the limit mg ^> T, where m# corresponds 
to the baryonic resonance mass, the Ansatz for the free energy density as a function 
of an imaginary chemical potential is 

~V¥ i Vf r = — VY 4 — = ^ - cos (^r)) > I 30 ) 

where f(T) = 4y J^ieBaryons (^r") 2 -^ (^r")- ^ e ^ us have a mean to measure the 
sum of resonances f{T). For example in the case of a 6 3 x 4 lattice, we find 
f(T ~ 0.9T C ) = 0.048(1) 8 . Our data can be well described by this Ansatz. This con- 
firms our expectation that the relevant degrees of freedom in the low-temperature 
phase are hadrons. The masses of these hadrons are much larger than the scale 
given by the temperature, since the free energy density changes only slightly when 
varying the imaginary chemical potential, thus m# ^> /xj ~ T c pa 160 MeV. 

For y ~ 1.1 we expect a cusp at fii = ±^y- (Z 3 -transitions) to develop in the 
thermodynamic limit, due to the first order phase transition. Indeed, it appears 
clearly as the volume increases, see Fig. El for a comparison of the histogram results 
(left) versus the reweighting approach 9 (right). 

We can try to describe these results by a generic Taylor series in as an Ansatz, 
which can be compared with a simple model at high temperature, the free gas of 



8 We thank D. Toublan |15j for estimating the sum of resonances for the four-flavour continuum 
theory. However, the result. f(T ~ 0.9T C ) rj 0.2, differs from our determination by about a factor 
4. It is unclear what is the main reason for this discrepancy, but the small, coarse lattice we use 
(a ~ 0.3 fm) certainly contributes an important part. 

9 Thc 8 3 x 4-data points are taken from the histogram method. 
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T/T„ - 1.1 



T/T„~1.1 





Figure 5: AF ^!^ as a function of ^- for ^ ~ 1.1. (Ze/£) The histogram method; 
[right) the reweighting method [3], supplemented by the histogram results for 8 3 x 4. 
A simple modification of the free gas expression describes all the data. As the volume 
increases, the data come close to the Stefan-Boltzmann limit (T — > oo) even though 
Y ~ 1.1 only. 



massless quarks. If we perform an analytic continuation from real to imaginary 
chemical potential, then the free energy density of this model is given by 

AF(T, = Nf_ (N\ 2 _ jV> (^V m 
VT A 2 \T/ 4tt 2 V T J ' 1 ' 

These simple expressions are valid in the continuum theory at very high temper- 
ature, where the coupling g(T) pa 0. On the lattice we expect finite size corrections 
(N s < oo) as well as cut-off corrections (T = Ref. [H)] has calculated the free 

energy of free fermions on a lattice having infinite spatial size (N s = oo) but finite 
temporal extent (N t = 4). Here, we also determine the corrections for finite spatial 
size N s = 4, 6, 8, 10 for the free massless fermion gas on the lattice. We set the gauge 
fields A^(x) = 0, ie. the gauge links to the identity, and solve for the free energy via 

AF^ e (T, fij) log Zf re %T, fn) log det M^(T, m) 



VT 4 VT 3 VT 3 

N f f^V n N f (ViV 



where C2 and C4 are fit coefficients. Table ^ summarises the results. 

The coefficients C2 and C4 approach their infinite volume expectation rather 
quickly. For the particular quark mass ^ = 0.2 which we consider, the difference 
from the massless limit is smaller than the (fitting) errors, and thus, results are not 
presented explicitly. Note that we have an additional column [Cq\: we have added 
the term C 6 (^) 6 to the Ansatz Eq. (p?2|) . The coefficient C 6 is very small and leaves 
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Lattice 


c 2 


c 4 


[C 6 ] 


4 3 x 4 


4.387(1) 


0.28(3) 


H 


6 3 x 4 


2.628(1) 


1.70(5) 


[0.0081(1)] 


8 3 x 4 


2.315(1) 


2.25(5) 


[0.0046(1)] 


10 3 x 4 


2.250(1) 


2.49(5) 


[0.0030(1)] 


oo 3 x 4 


2.25 


2.6 





Table 1: The prediction for the free energy density based on the free massless gas 
model in the continuum at high temperature suffers from finite size and cut-off 
effects. The correction terms C 2 and C4 help to quantify the systematics. The 
functional form in Eq. (j32|) is sufficient, since the contribution of the additional term 
(^) 6 is very small. 



C 2 and C 4 unchanged within the errors. 

In the end, we consider the volume-dependent lattice corrections C 2 and C4 and 
measure the deviation from this free gas model by two parameters b 2 (T) and 64 (T). 
The Ansatz to describe our results in Fig. 03 is thus 



AF(T, ///) 




T ~ 1.1T C 


k CO 


h(T) 


64 (T) (periodic) 


4 3 x 4 


0.29(1) 


32(4) 


7(1) 


6 3 x 4 


0.71(1) 


4.1(4) 


0.2(7) 


8 3 x 4 


0.90(2) 


3.6(8) 


1.4(4) 


SB limit (T -> 00) 


1 


1 


1 



Table 2: The coefficients of the free energy density expansion for ^- ~ 1.1 come 
close to their Stefan Boltzmann (T — > 00) value. There are two values for 64 (T): 
the first one is the result of the chi-square fit of Eq. (j3^j) : the second one ("periodic") 
makes use of a periodicised Ansatz, see Eq . (|34j) . The comparison of the two values 
gives some measure of the systematic error. 

We observe that the leading term approaches the Stefan Boltzmann limit rather 
fast upon increasing the volume, see Fig. El and Table El This is somewhat surprising 
since this coincidence with the Stefan Boltzmann law will occur only at T — > 00. 
Deviations at T ~ 1.1 T c should persist even in the thermodynamic limit, reflecting 
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the interactions of the quarks 10 . The reduction of 62 {T) from 1 is consistent with 
leading perturbative corrections ^7j. The value we obtain is consistent with that 
measured in Ref. ^H]- The simple prediction of the free massless quark gas model 
works better than expected. Thus, the relevant degrees of freedom at high tempera- 
ture are very light quarks, which is also visible in the strong dependency of the free 
energy density on the imaginary chemical potential, hence m q <C fii ~ T c « 160 
MeV. 

The coefficient 64 (T) in Table El suffers from systematic fitting errors. One source 
is the fitting range: our Ansatz Eq. (j33J) does not reflect the ^-periodicity of ^f, 
therefore we are allowed to fit small only. In this regime, the quartic term is 
subleading and hard to quantify. An estimate of the systematic fitting error can be 
obtained by varying the fitting range (not explicitly tabulated). Another source is 
the fitting Ansatz: we could add the next-order term (^r) 6 , which changes C 4 by a 
few percent, or periodicise the Ansatz by hand via 

A-F peT .(T, /ij) _ 1 , Z per (T, fij) 

VT 4 ~ VT 3 ° g Z per (T,0) 1 } 

with 

Z per (T^j) = jt e -^(^V^(^¥) 2 -MT)^(^¥) 4 ) . (35) 

k=— oo 

However, we must truncate (in practice, at k = ±1) the sum over all sectors to 
preserve convergence, because of the sign of the C4 contribution. In conclusion, we 
cannot determine 64 (T) accurately. Nevertheless, it is remarkable how well the free 
quark gas model describes our results. On an 8 3 x 4 lattice, the deviation is only 
about 10%. By using a canonical approach to simulate finite density QCD we 
can obtain more accurate results, to be presented in a follow-up publication. In 
particular, we can show that 6 4 (T = 1.1 T c ) = 2.97(3). 

7 Deconfinement Transition and Finite Size Ef- 
fects 

The phase transition is signaled by the peak in the susceptibility of the chiral con- 
densate (ipip) or in the specific heat. In Fig. El we show the former. 

10 It has been shown already that the free energy of the gluon sector deviates from the Stefan 
Boltzmann value by about 15% ^Hl even at ^- ~ 5 (and more for lower temperatures). It thus 
would be natural to observe deviations at finite temperature also in the quark sector. 
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Figure 6: Susceptibility of ipip versus (3 for all volumes in the grand canonical and 
canonical ensembles, left. Zqc(T^ = 0), right Zc(T,B = 0). Even for the 
smallest, 4 4 , lattice, differences are barely visible. 



On the 4 3 x 4 lattice, a slight shift in the pseudo-critical (3 C can be observed 
between the grand-canonical and the canonical results. It disappears for larger 
volumes. We observe the same behaviour for the specific heat. The small deviation 
is caused by contributions from B ^ sectors, which are present in the \x = 
ensemble. However they are suppressed by a factor ~ e~ ~ <^ \ } where is the 
mass of a baryon. In terms of the baryon density p, we recognise the exponential 
suppression in the volume since e" ~ = e~ p ~ . Thus, we verify once more 
that the zero chemical potential ensemble is equivalent to the zero baryon density 
ensemble in the thermodynamic limit. 

Note that the non-zero triality sectors have zero partition function and do not 
contribute. They do not affect observables studied in this section, which are insen- 
sitive to the centre symmetry. 

In quenched simulations, a Z 3 -symmetrisation of the Polyakov loop is sometimes 
enforced by hand, which is accompanied by reduced finite size effects [20] • Therefore, 
we might expect to observe a similar reduction in the canonical formalism as well 
compared to the grand canonical one. To compare the finite size effects in the two 
ensembles, we analyse the minimum of the Binder cumulant 



1 (O 4 ) 

C B (0) = 1 - (36) 

versus the inverse volume 1/V (see Fig. EJ). For both the plaquette and the chiral 
condensate, the thermodynamic (linear) extrapolation does not tend to | - indicative 
of a first order phase transition 11 , confirming the finding in the literature [22] for 

11 In the case of a second order transition, (O 4 ) is equal to (O 2 ) 2 up to finite size corrections |21j . 
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Figure 7: Binder cumulant minimum versus inverse volume for both ensembles 
(slightly shifted in the x-axis to enhance visibility), left: O =plaquette, right: 
O =chiral condensate. The thermodynamic extrapolation does not reach | (the 
upper left corner of the figure), indicating a first order transition. Finite-size effects, 
reflected in the 1 /V slope, are equivalent for both ensembles. 



our quark masses. However, for each volume, the measured cumulant values agree 
between the two ensembles within statistical errors, indicating equivalent finite size 
effects. 



8 Conclusions 

For all densities, volumes, and (finite) temperatures, the Polyakov loop expectation 
value is non-zero in the grand canonical ensemble Eq.(JJ), and zero in the equiva- 
lent canonical ensemble Eq. ([%]). This Polyakov loop paradox has to be considered 
an artifact of keeping, in the grand canonical ensemble, sectors with quark num- 
bers not multiple of three. These canonical sectors, the so-called non-zero triality 
sectors, have zero partition function. Thus, the non-vanishing expectation value 
(Pol) z Ga (T s n) i n the common grand canonical formulation of QCD at finite tempera- 
ture and density, Eq.([Tj). is irrelevant for thermodynamic properties. The physically 
meaningful Polyakov loop correlator (Pol(0)Pol(xY) behaves in the same way in 
both ensembles. 

Because of quantum and thermal fluctuations, (Pol(0)Pol(xy) tends to a non- 
zero value when |x| — > oo. On the other hand, (Pol) = in the canonical ensemble. 
Thus, the clustering property is violated, which shows that the center symmetry is 

Thus. Cb(0) — ► f in the thermodynamic limit. In the case of a first order transition, the double 
peak structure of the distribution of the measurements causes a non-trivial value of the Binder 
cumulant. 
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spontaneously broken in the canonical ensemble, rather than explicitly broken by 
the fermion determinant as in the usual grand-canonical ensemble. 

Furthermore, an explicitly centre-symmetric grand canonical partition function 
Z GC (T,fi), Eq.flSJ), can be constructed from the canonical partition functions, where 
the contributions of non-zero triality states are projected out. This partition function 
will give identical expectation values to the usual Zqc{T, fi), apart from a vanishing 
expectation value for the Polyakov loop. Therefore, the non-zero triality states can 
be included or excluded: the thermodynamic properties of the theory are unchanged. 

We have shown this explicitly by comparing the grand canonical ensemble at 
fi = with the canonical ensemble B — 0. Numerically, we have established the 
equivalence of Zqc{^ — 0) and Zc{B = 0) in the thermodynamic limit by measuring 
the free energy density as a function of fii, using the histogram of the imaginary 
chemical potential distribution 12 . For all temperatures, we observe a minimum at 

= 0. At low temperature, the free energy density of the confined phase can 
be rather well described by the hadron resonance gas, see Eq. (JHU|) . We thus have 
a simple way to determine the sum of resonances f(T). At high temperature, a 
slightly modified free gas Ansatz, see Eq. (jHH|) . allows to account for all data points 
in the quark-gluon plasma phase. We determine the finite-size and cut-off correction 
terms C2 and C4 by calculating the free energy of the free fermion gas on the lattice 
and find agreement with the literature for V — ■> 00. By construction, the interaction 
coefficients b 2 (T) and 64 (T) tend to 1 for high enough temperatures, reproducing the 
Stefan Boltzmann law. Just above T c , the deviation from 1 in the leading coefficient 
is about 30% on a 6 3 x 4 lattice; on an 8 3 x 4 lattice, this deviation is about 10% only. 
This near-agreement with a non-interacting gas is unexpected at such comparatively 
low temperatures. 

The approach to the thermodynamic limit is very similar in the canonical (B = 0) 
and grand canonical (jj, = 0) ensembles. The susceptibility of the chiral condensate, 
or the specific heat, indicate the same pseudo-critical temperature already on small 
volumes. A small shift, caused by contributions from non-zero baryon sectors, is 
visible only on the 4 3 x 4 lattice. 

The zero-density canonical formulation requires a centre-symmetric simulation 
of QCD, which can be achieved very simply with negligible computer overhead, by 
adding to the standard algorithm a single degree of freedom fij updated by Metropo- 
lis. 

We hope to have fully clarified the (un) importance of non-zero triality states, and 
12 Rcmcmber that our numerical approach treats [ii as a dynamical degree of freedom. 
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thus, to have put to rest long-standing speculations. Further connections between 
the grand-canonical and the canonical formalisms in the context of non-zero chemical 
potential/density will be the subject of a forthcoming paper. 
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Appendix: Stochastic Estimator 

A ratio of determinants can be estimated using a single Gaussian complex vector: 



det^(ff(/4) +m) _ det*' M(/^) _ / dfidfa M f W 
det N f(p (fij) +m) ~ det N f M(//j) " 



_ fdrfd V 17(0,^,^)1 e -^M"f'\» I )M-»f'\»> I )M~»f'\»>w N i /2 ^)v 
J drfdrj | J(0, //,///) | e-^" 

(38) 

where we have substituted <p = M Nf / 2 ((ii)r]. Note, that in the above notation, the 
Jacobian \J(<f), Vi^i)] is det M N f (fij), which is independent of rj and cancels out in 
the ratio: 

det N f (lp{Lii) + m)~ J drfdrj e~^ 2 ^ 

(■)r) tells us that r] has to be sampled with the distribution J drfdr/ e~^ 2 . 
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